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ACCELERATING STOCHASTIC COLLOCATION METHODS FOR PARTIAL 
DIFFERENTIAL EQUATIONS WITH RANDOM INPUT DATA * * * § 

D. GALINDOt, P. JANTSCHt, C. G. WEBSTER?, AND G. ZHANG’ 

Abstract. This work proposes and analyzes a generalized acceleration technique for decreasing the computational 
complexity of using stochastic collocation (SC) methods to solve partial differential equations (PDEs) with random 
input data. The SC approaches considered in this effort consist of sequentially constructed multi-dimensional La¬ 
grange interpolant in the random parametric domain, formulated by collocating on a set of points so that the resulting 
approximation is defined in a hierarchical sequence of polynomial spaces of increasing fidelity. Our acceleration ap¬ 
proach exploits the construction of the SC interpolant to accelerate the underlying ensemble of deterministic solutions. 
Specifically, we predict the solution of the parametrized PDE at each collocation point on the current level of the SC 
approximation by evaluating each sample with a previously assembled lower fidelity interpolant, and then use such 
predictions to provide deterministic (linear or nonlinear) iterative solvers with improved initial approximations. As a 
concrete example, we develop our approach in the context of SC approaches that employ sparse tensor products of 
globally defined Lagrange polynomials on nested one-dimensional Clenshaw-Curtis abscissas. This work also provides 
a rigorous computational complexity analysis of the resulting fully discrete sparse grid SC approximation, with and 
without acceleration, which demonstrates the effectiveness of our proposed methodology in reducing the total num¬ 
ber of iterations of a conjugate gradient solution of the finite element systems at each collocation point. Numerical 
examples include both linear and nonlinear parametrized PDEs, which are used to illustrate the theoretical results 
and the improved efficiency of this technique compared with several others. 

Key words, stochastic and parametric PDEs, stochastic collocation, high-dimensional approximation, un¬ 
certainty quantification, sparse grids, multivariate polynomial approximation, iterative solvers, conjugate gradient 
method 
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1. Introduction. Modern approaches for predicting the behavior of physical and engineering 
problems, and assessing risk and informing decision making in manufacturing, economic forecasting, 
public policy, and human welfare, rely on mathematical modeling followed by computer simulation. 
Such predictions are obtained by constructing models whose solutions describe the phenomenon of 
interest, and then using computational methods to approximate the outputs of the models. Thus, 
the solution of a mathematical model can be viewed as a mapping from available input information 
onto a desired output of interest; predictions obtained through computational simulations are merely 
approximations of the images of the inputs, that is, of the output of interest. There are several causes 
for possible discrepancies between observations and approximate solutions obtained via computer 
simulations. The mathematical model may not, and usually does not, provide a totally faithful 
description of the phenomenon being modeled. Additionally, when an application is considered, 
the mathematical models need to be provided with input data, such as coefficients, forcing terms, 
initial and boundary conditions, geometry, etc. This input data may be affected by a large amount 
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of uncertainty due to intrinsic variability or the difficulty in accurately characterizing the physical 
system. 

Such uncertainties can be included in the mathematical model by adopting a probabilistic set¬ 
ting, provided enough information is available for a complete statistical characterization of the 
physical system. In this effort we assume our mathematical model is described by a partial differ¬ 
ential equation (PDE) and the random input data are modeled as finite dimensional random fields, 
parameterized by a vector y = (j/i, • • • , yjsi) of dimension TV, consisting of uncorrelated real-valued 
random variables. Therefore, the goal of the mathematical and computational analysis becomes the 
approximation of the solution map y i—>■ u(y), or statistical moments (mean, variance, covariance, 
etc.) of the solution or some quantity of interest (Qol) of the system, given the probability dis¬ 
tribution of the input random data. A major challenge associated with developing approximation 
techniques for such problems involves alleviating the curse of dimensionality, by which the compu¬ 
tational complexity of any naive polynomial approach will grow exponentially with the dimension 
N of the parametric domain. 

Monte Carlo (MC) methods (see, e.g., 17|) are the most popular approaches for approximating 
high-dimensional integrals, based on independent realizations u{yk), k = 1,... ,M, of the parame¬ 
terized PDE; approximations of the expectation or other Qols are obtained by averaging over the 
corresponding realizations of that quantity. The resulting numerical error is proportional to 
thus achieving convergence rates independent of dimension N, but requiring a very large number 
of samples to achieve reasonably small errors. Other ensemble-based methods, including quasi-MC 
(QMC) and important sampling (see [24, 2^, 3^ and the references therein), have been devised to 
produce increase convergence rates, e.g., proportional to log(M)’’(^\ however, the function 
r{N) > 0 increases with dimension N. Moreover, since both MC and QMC are quadrature tech¬ 
niques for Qols, neither have the ability to simultaneously approximate the solution map y i—>■ u{y), 
required by a large class of applications. 

In the last decade, two global polynomial approaches have been proposed that often feature much 
faster convergence rates: intrusive stochastic Galerkin (SC) methods, constructed from pre-defined 
orthogonal polynomials [3,13, or best M-term and quasi-optimal approaches S i B El , and 
non-intrusive stochastic collocation (SC) methods, constructed from (sparse) Lagrange interpolating 


polynomials [H, [3) 311, or discrete projections 27 


. These methods exploit the underlying 
regularity of the PDE solution map u{y) with respect to the parameters y, evident in a wide class 
of high-dimensional applications, to construct an approximate solution, and differ only in the choice 
of basis. 

For both SG and SC approaches, the overall computational cost grows rapidly with increasing 
dimension. A recent development for alleviating such complexity and accelerating the convergence 
of parameterized PDE solutions is to utilize multilevel methods (see e.g., multilevel Monte Carlo 
(MLMC) methods @,[1 El 013 and the multilevel stochastic collocation (MLSC) approach [4lj|^- 
The main ingredient to multilevel methods is the exploitation of a hierarchical sequence of spatial 
approximations to the underlying PDE, which are then combined with discretizations in parameter 
space in such a way as to minimize the overall computational cost. The approximation of the solution 
u on the finest mesh is represented by the approximation on the coarsest mesh plus a sequence of 
“correction” terms. The resulting decrease in complexity with the use of multilevel methods results 
from the fact that the dominant behavior of the solution u can be captured with cheap simulations 
on coarse meshes, so that the number of expensive simulations computed on fine meshes can be 
considerably reduced. 

Nonetheless, the dominant cost in applying any uncertainty quantification (UQ) approach lies 
in the solution of the underlying parametrized linear/nonlinear PDFs, for a given value of the 
random inputs. Such solutions are often computed using iterative solvers, e.g., conjugate gradient 
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(CG) methods for symmetric positive-definite linear systems, generalized minimal residual method 


(GMRES) for non-symmetric linear systems [37|, and fixed-point iteration methods[36| for nonlinear 
PDEs. However, many high-fidelity, multi-physics models can exhaust the resources of the largest 
machines with a single instantiation and, as such, are not practical for even the most advanced UQ 
techniques. As such, several methods for improving the performance of iterative solvers have been 
proposed; especially preconditioner and subspace methods for iterative Krylov solvers. A strategy 
that utilizes shared search directions for solving a collection of linear systems based on the GG 


method is proposed in [8|. In [33| a technique called Krylov recycling was introduced to solve sets of 
linear systems sequentially, based on ideas adapted from restarted and truncated GMRES (see [s^ 
and the references therein). This approach was later applied to the linear systems that arise from 
SG approximations that use the so-called doubly orthogonal bases to solve stochastic paramterized 
PDEs 2^ . In addition, several preconditioners have been developed that improve the performance of 
solving the large linear systems resulting from SG approximations that employ standard orthogonal 

polynomials a Safi. 

On the other hand, when a general linear solver is employed to solve the underling SG or SG 
approximation, it is straightforward to see that improved initial approximations can significantly 
reduce the number of iterations required to reach a prescribed accuracy. A sequential orthogonal 
expansion is utilized in Bin such that a low resolution solution provides an initial guess for the 
solution of the system with an enriched basis. However, at each step, all the expansion coefficients 
must be explicitly recomputed, resulting in increased costs. Similarly, in [2lj an extension of a 


mean-based preconditioner is applied to each linear system coming from a sequential SG approach, 
wherein the solution of the j-th system is given as the initial vector for the (j -I- I)-th system. This 
approach, as well as the Krylov recycling method, impose an ordering of the linear systems that 
appear in the SG approximation. Gonsequently, new approaches are needed to amortize the cost 
of expensive simulations by reusing both deterministic and stochastic information across multiple 
ensembles of solutions. 

In this work, we propose to improve the computational efficiency of non-intrusive approxima¬ 
tions, by focusing on SG approaches that sequentially construct a multi-dimensional Lagrange inter- 
polant in a hierarchical sequence of polynomial spaces of increasing fidelity. As opposed to multilevel 
methods that reduce the overall computational burden by taking advantage of a hierarchical spatial 
approximation, our approach exploits the structure of the SG interpolant to accelerate the under¬ 
lying ensemble of deterministic solutions. Specifically, we predict the solution of the parametrized 
PDE at each collocation point on the current level of the SG approximation by evaluating each 
sample with a previously assembled lower fidelity interpolant, and then use such predictions to pro¬ 
vide deterministic (linear or nonlinear) iterative solvers with improved initial approximations. As a 
particular application, we pose this acceleration technique in the context of hierarchical SG methods 
that employ sparse tensor products of globally defined Lagrange polynomials Blill , on nested 
one-dimensional Glenshaw-Gurtis abscissas. However, the same idea can be extended to other non- 
intrusive collocation approaches including orthogonal polynomials , as well as piecewise local and 
wavelet polynomials expansions UM- 

The sparse grid SG approximation considered in this work produces a sequence of interpolants, 
where a new set of collocation points is added on each level in order to increase the accuracy of the 
interpolant. For each newly added collocation point on the current level, we predict the solution of 
the underlying deterministic PDE using the most up to date sparse grid interpolant available; the 
previous level’s interpolant. We then use the prediction as the starting point of the iterative solver. 
The uniform convergence of the sparse grid interpolant to the true solution results in an increasingly 
accurate initial guess as the level increases, so that the overall complexity of the SG method can 
be dramatically reduced. We apply our novel approach in the context of solving both linear and 
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nonlinear stochastic PDEs, wherein, we assume that the parameterized systems are solved by some 
existing linear or nonlinear iterative method. Furthermore, in the linear case, this technique can 
also be used to efficiently generate improved preconditioners for linear systems associated to the 
collocation points on higher levels, which further accelerates the convergence rate of the underlying 
solver. 

The outline of this paper is as follows: We begin by describing the class of parameterized 
linear and nonlinear stochastic PDEs under consideration in ^ In ^we describe our acceleration 
technique in the context of general stochastic collocation methods, defined on a hierarchical sequence 
of polynomial spaces, for approximating both linear and nonlinear stochastic elliptic PDEs using 
nonlinear iterative solvers. In ^ we briefly recall the sparse grid SC method, where the sparse 
grid interpolant is constructed with the use of nested one-dimensional Clenshaw-Curtis abscissas. 
The theoretical convergence rates, with respect to the level of the interpolant and the degrees of 
freedom are shown in 114.II In 1 14.21 we provide a rigorous computational complexity analysis of 
the resulting fully discrete sparse grid SC approximation, with and without acceleration, used to 
demonstrate the effectiveness of our proposed methodology in reducing the total number of iterations 
of a conjugate gradient solution of the finite element systems at each collocation point. Finally, in 
^ we provide several numerical examples, including both moderately large-dimensional linear and 
nonlinear parametrized PDEs, which are used to illustrate the theoretical results and the improved 
efficiency of this technique compared with several others. 

2. Problem setting. Let D C = 1,2,3, be a bounded domain and let (D,J^,P) denote 
a complete probability space with sample space D, cr-algebra J- = 2^, and probability measure 
P : — >■ [0,1]. Define £ as a differential operator that depends on a coefficient a{x,u!) with x G D 

and uj G fl. Analogously, the forcing term / = f(x,Lu) can be assumed to be a random field as well. 
In general, a and / belong to different probability spaces but, for economy of notation, we simply 
denote the stochastic dependences in the same probability space. Consider the stochastic boundary 
value problem. Find a random function u : D x D —>■ K. such that, P-a.e. in D, the following equations 
hold: 


I ^a){u) = f in D, 

I u = g on dD, 

where g is a suitable boundary condition. We denote by W{D) a Banach space and assume the 
underlying random input data are chosen so that the corresponding stochastic system (EH) is well- 
posed and has a unique solution u{x,u!) G Lp(fl; W(-D)), the function space given by 


L^(D;W(D)) 


u : D X fl 


u is strongly measurable and 


< +C)0 


In this setting, the approximation space consists of Banach-space valued functions that have finite 
g-th order moments. Two example problems posed in this setting are given as follows. 

Example 2.1. (Linear elliptic problem). Find a random field m : D x D — > R such that P-a.e. 


J — V • (a(a;, a;)Vu(a;, w)) 
1 u{x,uj) 


= 

= 0 


in D X fl, 
on dD X D, 


( 2 . 2 ) 
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where V denotes the gradient operator with respect to the spatial variable x € D. The well-posedness 
of (12.21) is guaranteed in Lp{fl; Hq{D)) with a(x,uj) uniformly elliptic, i.e., 


P(c^ G Tl O-inin ^ a{x,Uj') ^ C^max Vx € — 1 with Ctmin; ^max 

and f{x,uj) square integrable, i.e., 

[ E[/^]c?x := f ( p{x,ui) dF{ui)dx <+CO. 

Jd Jd Jq 


€ (O,oo), 


(2.3) 


Example 2.2. (Nonlinear elliptic problem). For k G N, find a random field u : D x Tl ^ 
such that P-a.e. 


—'S/-{a{x,uj)'S/u{x,uj))+u{x,u!)\u{x,uj)\^ = f{x,Lo) in D, 

u{x, w) = 0 on dD. 


(2.4) 


The well-posedness of (1^ is guaranteed in Lp (fl; W(D)) with a, f as in Examvle \2.1\ and W(D) = 

In many applications, the source of randomness can be approximated with only a finite number 
of uncorrelated, or even independent, random variables. For instance, the random input data a and 
/ in p.lll may have a piecewise representation, or in other applications may have spatial variation 
that can be modeled as a correlated random field, making them amenable to approximation by a 
Karhunen-Loeve (KL) expansion (^ . In practice, one has to truncate such expansions according to 
the desired accuracy of the simulation. As such, we make the following assumption regarding the 
random input data a and / (cf [2^|3l|). 

Assumption 2.1. (Independence and finite dimensional noise). The random fields a(x,uj) and 
f{x,Lu) have the form: 

a{x,ui) = a{x,y{uj)) and f(x,Lo)=f{x,y(Lo)) on D x Q, 

where yfjj) = [yi(a;),... ,y]sr{ui)] : —>■ is a vector of independent and uncorrelated real-valued 

random variables. 

We note that Assumption 12.11 and the Doob-Dynkin lemma guarantee that a{x,y{ui)) and 
f{x,y{uj)) are Borel-measurable functions of the random vector y : ft ^ R^. In our setting, we 
denote by r„ = 2/„(f2) C M the image of the random variable ?/„, and set F = nr!,=i where 
N G N+. If the distribution measure of y{uj) is absolutely continuous with respect to Lebesgue 
measure, then there exists a joint probability density function of y(w) denoted by 

N 

g{y):T^R+, with g{y) = Y[ QuiVn) & (T)- 

n—1 

Therefore, based on Assumption l2.11 the probability space (II, J", P) is mapped to (r,S(F), g[y)dy), 
where S(r) is the Borel cr-algebra on F and g(y)dy is a probability measure on S(r). By assuming 
the solution u of mi is (T-measurable with respect to a and /, the Doob-Dynkin lemma guarantees 
that u{x,uj) can also be characterized by the same random vector y, i.e.. 


u[x,u}) = u{x,yi{uj),...,yN{aj)) G L«(F;IF(D)), 
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where L|(r;Vh(D)) is defined by 


Lg(r; VF(Z))) = |M:Dxr —u strongly measurable and J ||'u||^(£)) Q{y)dy < oo 
Note that the above integral will be replaced by the essential supremum when q = oo: 


L°°{r;W{D)) 


|it: I? X r —>• 


u strongly measurable and esssup^ \\u{y)\\^^f^jy^ 


< oo 


2.1. Weak formulation. In what follows, we treat the solution to m as a parameterized 
function u{x.Ji) of the A^-dimensional random variables y G T C . This leads to a Galerkin weak 
formulation [23j of the PDE in (EH), with respect to both physical and parameter space, i.e., seek 
u G T|(r; W(D)) such that 



;y)T^{v) gdxdy 


fvgdxdy, Vu e L«(r; W(£))), 


r JD 


where T^ji' G Ai U A 2 are linear operators independent of y, while the operators Si, are linear for 
ly G All s-nd nonlinear for v G A 2 . Moreover, since the solution u can be viewed as a mapping 
M : r —^ W{D), for convenience we may omit the dependence on x G D and write u(y) to emphasize 
the dependence of u on y. As such, we may also write the problem (j2.1|) in the alternative weak 
form 

f { ^ 'S'!^('«(y);y)Th(u) ) dx = / f{y)vdx, Wv gW{D), g-a.e. inT. (2.5) 

Vt-GAiUAa / 


Therefore, the stochastic boundary-value problem EH has been converted into a deterministic 
parametric problem (12.51) . The acceleration technique proposed in (Eland the sparse-grid SC method 
discussed in Elwill be based on the solution of the weak form (12. 5|) above. 

3. Accelerating stochastic collocation methods. Our acceleration scheme will be proposed 
in the context of both linear and nonlinear elliptic PDEs. A general SC approach requires the semi¬ 
discrete solution Uh{-,y) G Wh{D) C W{D) at a set of collocation points {yn j}^i C P, given 
by 


Mh 

Uh{x,yL,j) ='^CL,j,i(Pz{x), j = l,...,Mi. (3.1) 

Here is a predefined finite element basis of WhiD), and for j = 1,..., M^, the coefficient 

vector clj := {clj,!, ■ ■ ■ ,CL,j,Mh)^ is the solution of the following system of equations: 


Sy{ifi-,yL,j) T„{Lpi,)dx 

i=l i^gAi 

f /Mu \ 

= fiyL,j)^^'- ^ s„ l'^CL,j,i(pi-,yL,j]T„{(pi,) dXi i' = l,...,Mhi 

i'GAs \i=l / 


(3.2) 
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with Si, and T„ defined as above. Note that is equivalent to (12.51) with the nonlinear oper¬ 
ators subtracted on the right hand side. When A 2 = 0, the PDE is linear, and a standard FEM 
discretization leads to a linear system of equations. 

For L G N+, we denote by II an interpolation operator that utilizes Ml collocation points, 
defined hy 'Hl = More generally, assume that we have a family of interpolation operators 

which for each L G N+ approximates the solution Uh{x, •) in polynomial spaces 

iPi(r) c ... c VUT) c VL+iir) c ... c Ll(r), 


of increasing fidelity, defined on sets of sample points TLl C E. Assume further that the fully discrete 
solution Uh^L G Wh{D) ® PL(r) has Lagrange interpolating form 


Ml /Mh \ 

Uh,L{x,y) := lL[uh]{x,y) i'^CL,j,z(Pi{x)\ 'i’L,j{y), (3.3) 

j=i \i=i / 

where is a basis for The approximation (13.3|) can be constructed by solving for 

Uh{x,yL,j) independently at each sample point y^j G Hl- In (31 we construct a specific example of 
an interpolation scheme satisfying (IXX) . namely global sparse grid collocation. 

For each L G N, the bulk of the computational cost in using (13.31) goes into solving the Ml 
systems of equations (13.21) corresponding to each collocation point yL,j, j = 1, ■ ■ ■, Ml- Since the 
systems are independent and deterministic, they can be solved separately using existing FEM solvers, 
providing a straightforward path to parallelization compared to intrusive methods such as stochastic 
Galerkin methods. In this work, we consider iterative solvers for the system in (13.21) . and propose an 
acceleration scheme to reduce the total number of iterations necessary to the collection of systems 
over the set of sample parameters. 

Denoting by ut the output of the selected iterative solver for the system (13.21) . for yLy G Hl 
the semi-discrete solution Uh{x,yL,j) is approximated by 


Mh Mh 

Uh{x,yL,j) = '^CLj,iPi{x) Ki dh{x,yL,j) = Pi{x), 

where we define clj = (cl.j.i, ■ • ■ ,CL,j,Mh)^^ therefore the final SC approximation is given by 
a perturbation of (13.31) . i.e.. 


Ml /Mh \ 

Uh,L{x,y) ■■= ^ i'^CL,j,iipi{x) I (3.4) 

j=l Vi=l / 

We observe that the performance of the underlying iterative solver can be improved by proposing 
a good initial guess, denoted or constructing an effective preconditioner to reduce the condition 
number of the system. Here, we propose our approach for improving initial deterministic approxima¬ 
tions, remarking that the same idea can be also utilized to construct preconditioners. To start the 
iterative solver for the system in (13.21) . it is common to use a zero initial guess, i.e., = (0,..., 0)^. 

However, we can predict the solution at level L using lower level approximations to construct im¬ 
proved initial solutions c^l\- Assume that we first obtain Uh^L-iix,y) by collocating solutions to 
(IQ) over Hl- 1 - Then at level A, for each new point yL,j GHl \ Hl-i, the initial guess can 
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be given by interpolating the solutions from level L — i.e., 

<^^,3 = {^h,L-l{Xi,yL,]),■■■,Uh,L-l{XMh,yL,j)^ = ^ L-l,j'{yL,j)- (3.5) 

J' = l 

For a convergent interpolation scheme, we expect the necessary number of iterations to compute 
CLj to become smaller as the level L increases to an overall maximum level, denoted Fmax- As such, 
the construction of the desired solution Uh,LmB.^ is accelerated through the intermediate solutions 
{^h,L}^=‘T~^■ Note that this approach reduces computational cost by improving initial guesses, 
but does not depend on the specific solver used. Thus, our scheme may be combined with other 
techniques for accelerating convergence, such as faster nonlinear solvers or better preconditioners. 
When the underlying PDF is nonlinear with respect to u, iterative solvers are commonly used for 
the solution of (j3.2L In Algorithm 1, we outline the acceleration procedure described above, using 
a general nonlinear iterative method for the solution of (13.2L 


Algorithm 1: The accelerated SC algorithm 
Goal: Compute Uh,L^^^{x,y) := 

1 : Define Mq = 1 and co,i = (0,..., 0)^ 

2; for L = 1, . . . ,Lmax do 

3; for yL,j &Ul \ '^/) do 

4: Compute the initial guess according to (1331); 

6: Initialize: k = 1 

7: repeat 

8: Compute residual = (r^L^,!^ ■ ■ ■. ' 

9: for z = I,..., M/j do 


10 : 

11 : 

12 : 

13: 

14: 

15: 

16: 


Jk) _ 


Id f (yL,j) Ti - 4y.i' Ti' (x), yL,]^ dx 


end for 

Update the solution; ■ 
k = k + 1 


1/GA1UA2 


. ^(fc+i) _ (fc) 




(k) N 

LJ) 


until <r 


end for 


17: end for 


The efficiency of the proposed algorithm will depend crucially on the number of times the 
iterative solver is utilized, i.e., how many sample points are in the set A.TLl = TLl \ for 

each level L. In fact, if the sample points are not nested, it could be the case that /S.TLl ='Hl, and 
the algorithm may be very inefficient. Hence, in the following sections we will assume: 

Assumption 3.1. Assume that the multidimensional point sets 'Hl, A = 1,..., Umax xre nested, 





Hi c c ... c c r. 

Then ATLl = T~Ll \ TLl-i, and we can construct the intermediate solutions using a 

subset of the information needed to approximate M/i,Lmax- 

In ^ we construct an interpolant using a point set which satisfies Assumption 13.11 Next, we 
give several examples of Algorithm 1, using iterative solvers for both nonlinear and linear elliptic 
PDEs. 

Example 3.1. Consider the weak form of the nonlinear elliptic PDE in Examvle \2.2l letting 
Si{v, y) = a{x, y)Vw, Ti(v) = Vu, S 2 {v, y) = v{x, y)\v{x, y)!®, and T 2 {v) = v (note that Ai = {1}, 
A 2 = {2} ). When using the fixed point iterative method in Algorithm 1, for the update step we define 




L,j’ 


r-gj.) = 




where the matrix = A{yLj), j = 1,... ,Ml is defined by 


[ALj]i,i'= a{yL,j)V(pi>W(pidx, for i,i'= 1,..., Mh- 
J D 


(3.6) 


With yL,j) = i ^^*■5 update is eguivalent to solving the following linear system 


f a(yLj)v4^2^^ Vw dx = f \f{yL,j)-ul^^^{yL,j)\ul^iiyLjW]vdx Vv eWhiD), 
JD Jd ^ 


to update to at the (k-\-l)-th iteration. Note that each iteration of the solver in Algorithm 

1 requires the solution of this linear system, which is not accelerated by our algorithm. 

Example 3.2. As a special case of the example above, consider the weak form of the linear 
elliptic problem in Examvle \2.1\ with Ai = {1}, A 2 = 0, Si{v;y) = oN/v and Ti{v) = Vu in (13.211 . 
Due to the linearity, at each collocation point the solution Uh{x,yL,j) = can be 

approximated by solving the following linear system 


— fL,ji 


(3.7) 


with Alj = A{yL,j), j = as in ((3?6)) . and {fL,j)i = Jd f{x,yL,j)Ti{x)dx for i = 

1,..., Mh. Under our assumptions on the coefficient a, the linear system (IS31) is symmetric positive 
definite, and we can use the CG method to find its solution. For k G N"*", by recursively defining 


(k) (k) 

plj = n!, 


(/c')T A (k) 

Plj j (fe') 

■ ^ 

k'<kPLj 


we get the update function 




( 1 ) 


(fe) 

Cr 

L,3 


) = 


(fc)T 


PlY ^l,jP^l:3 


,(fc) 

Ef_(k) 
(k) ALj- 
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Recall the following well-known error estimate for CG: 


(fc) 

clj - Clj 


< 2 


\y/^L,j + 1 


( 0 ) 

CLj - eft 


Al,' 


(3.8) 


where = K{yL,j) denotes the condition number of cf ■ is the vector of initial guess and 

is the output of the k-th iteration of the CG solver. As opposed to Examvle \S.l[ for this example 
Algorithm 1 accelerates the solution of the linear system (13.71) . 

To evaluate the efficiency of the accelerated SC method, we define cost metrics for the construc¬ 
tion of standard and accelerated SC approximations. In general, the computational cost in floating 
point operations (flops) is the total number iterations to solve (13.21) summed over each of the sample 
points—denoted by ifzero and ifacc for the standard and accelerated SC methods, respectively— 
multiplied by the cost of performing one iteration, denoted Citer- Let Cint be the additional cost of 
interpolation incurred by using the accelerated initial vectors (1531) . Then, we define 


CzorO - Cifor K7.1 


for the standard SC approach, and 


Cacc — Gii 


■ -l- Cj 


int ? 


(3.9) 


(3.10) 


for the accelerated SC approximation, respectively. 

In Example 13.21 the discretization of the linear PDE leads to Ml sparse systems of equations of 
size Mh X Mh- When solving these systems with a CG solver, and i^acc are the sum of solver 

iterations contributed from each sample system. In this case, the cost of one iteration is just the 
cost of one matrix vector product, i.e., Citer = CdMu, where Co depends on the domain D and the 
type of hnite element basis. 

Remark 3.1. (Relationship to multilevel methods). Multilevel methods reduce the complexity 
of stochastic sampling methods by balancing errors and computational cost across a sequence of 
stochastic and spatial approximations. Let Uh^ S 14, k = 0,...,K, be a sequence of semi-discrete 
approximations built in nested spaces, i.e., Vq C ... C Vk- Multilevel methods are based on the 
following identity: 


K 

'^hK — 'y ~ 'a/tfc_i)- 

fc =0 

Letting Qr^-k-: k = 0,..., K, denote the chosen method of stochastic approximation, a general mul¬ 
tilevel method might be written as 


^{ML) ^ J. 

fc =0 

The main idea is that highly resolved, expensive stochastic approximations, e.g., Qlkj combination 
with coarse deterministic approximations, that is, Uhg, and vice versa. In a similar way, collocation 
with nested grid points provides a natural multilevel hierarchy which we use in our method to ac¬ 
celerate each PDE solve (liT31) . A combination of these methods could involve using our algorithm 
to accelerate the construction of the operators QhK-ki reusing information from level to 

level, thus improving further the performance of SG methods. 
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Remark 3.2. (Interpolation costs). Note that many adaptive interpolation schemes already 
require evaluation of the intermediate interpolation operators as in (1^ . e.g., to compute residual 
error estimators. Thus, these methods will incur the interpolation cost Cint even in the case of zero 
initial vectors. Furthermore, for most nonlinear problems the deterministic solver is expensive, thus 
reducing the number of iterations is the most important element in reducing the cost. In each of 
these settings, we can define the cost metrics to simply be K^em cmd RTacc- 

Remark 3.3. (Hierarchical preconditioner construction). When solving linear systems using 
iterative methods, convergence properties can be improved by considering the condition number of 
the system. As with initial vectors, an interpolation algorithm can be used to construct good, cheap 
preconditioners. We consider preconditioner algorithms where an explicit preconditioner matrix, or 
its inverse, is constructed. In this case, for some low collocation level Lpc, we construct a strong 
preconditioner, Rlpcj' •= each individual iterative solver, j = l,...,Mipp. Then, 

these lower level preconditioners are interpolated for the subsequent levels. More specifically, for 
L > Lpc, and yL,j & 'Hl \ Hlpc^ '>^e use the preconditioner 

Mppc 

Pl,] ■■= P{yL,j) = ^ PlpcT ^LpcfiVL.j)- (3.11) 

j'=l 

Numerical illustrations of this approach are given in 

4. Applications to sparse grid stochastic collocation. In this section, we provide a specific 
example of an interpolation scheme satisfying the assumptions described in lj31 i.e., a generalized 
sparse grid SC approach for a fixed level L. In what follows, we briefly review the construction of 
sparse grid interpolants, and rigorously analyze the approximation errors and the complexities of 
both the standard and accelerated SC approaches, in order to demonstrate the improved efficiency 
of the proposed acceleration technique when applied to iterative linear solvers. 

The fully discrete SC approximation is built by polynomial interpolation of the semi-discrete 
solution Uh{x, y) on an appropriate set of collocation points in T. In our setting, such an interpolation 
scheme is based on a sparse tensor products of one-dimensional Lagrange interpolating polynomials 
with global support. Specifically, in the one-dimensional case, iV = 1, we introduce a sequence of 
Lagrange interpolation operators : (^^(r) —>• with Pj„(z)-i(r) the space of degree 

m{l) — 1 polynomials over L. Given a general function v € C'*^(r), these operators are defined by 


m(l) 

Here I G N represents the resolution level of the operator, m(/) € N+ denotes the number of 
interpolation points on level I, ifl{y) = 1 and for 1 > 1, 

m(i) ^ ; 

ti yj - y' 


are the global Lagrange polynomials of degree m{l) — l associated with the point set'd* = {y {,..., 

To satisfy Assumption 13.II we need nestedness of the one-dimensional sets, i.e., C d’", which is 
determined by the choice of interpolation points and the definition of m(l). In addition, we remark 
that similar constructions for ean be built based on wavelets 22| or other locally supported 
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polynomial functions 23 1 . 

In the multi-dimensional case, 
the difference operator 


i.e., > 1 , using the convention that = 0 , we introduce 


(g) ••• 0 


N 

0 ^ 


n=l 


(4.1) 


and define the multi-index l=(/i,...,/Ar)£N(^. The desired approximation is defined by a linear 
combination of tensor-product operators (14.111 over a set of multi-indices, determined by the condition 
( 7 ( 1 ) < L, for L G N+, and ( 7 ( 1 ) : N+ a strictly increasing function. For v G (^“(r) , we now 

define the generalized SC operator I™’® by 


0---0 H ( y ) 

s(i)<i 

= H (-l)l'l H(y), 

g(l)<LiG{0,l}'^ 


(4.2) 


where i = (ii,..., ijv) is a multi-index with G {0,1}, |i| = ii + ■ ■ ■ In, and L G N+ represents 
the approximation level. This approximation lives in the tensor product polynomial space given by 

= span Y[ 

where the multi-index set is defined as follows 



A 


m,g 

L 


|igN^ 


y(m^(l + 1 )) 



Here m(l) = ..., m{lN)), and m^l) := min{ry G N+ : m{w) > Z} is the left inverse of m (see 

i)- 

Specific choices for the one-dimensional growth rate m(l) and the function 17 ( 1 ) are needed to 
define the multi-index set A™’® and the corresponding polynomial space for the approximation. 

In this work, we construct the interpolant in (14.2p using the anisotropic Smolyak construction, i.e.. 


N 

m(l) = 1, m{l) = 2^“^ -I- 1 for / > 1 and y(l) = {In — 1), (4.3) 

, ^min 
n—1 


where cx = (ai,..., ajv) G K+ is a vector of weights reflecting the anisotropy of the system, i.e., 
the relative importance of each dimension, with Omin := min„ q;„ (see for more details). Our 
analysis does not depend strongly on this choice of m and y, and we could use other functions, e.g., 
m{l) = I and y(l) = max„ define the anisotropic tensor product approximation. 

When r is a bounded domain in a common choice is the Clenshaw-Curtis abcsissas 0 
given by the sets of extrema of Chebyshev polynomials including the end-point extrema. For a 
sample set of any size m{l) > 1 , the abscissas in the standard domain [— 1 , 1 ] are given by 


= <1 y' G [-1,1] 


y] = - cos 


• (J - 1 ) 

.(0 - 1 


for j = 1 , 


,m (0 


(4.4) 
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By taking y\ = Q and letting m{l) grow according to the rule in (14.3L one gets a sequence of 
nested sets d* C for I G N+. In addition, with ff(i) defined as in dUl), the resulting set of 
A^-dimensional abscissas is a Clenshaw-Curtis sparse grid. Other nested families of sparse grids can 
be constructed from, e.g., the Leja points [l^, Gauss-Patterson ^|, etc. 

Remark 4.1. (Specific Choice of m,g). For the remainder of the paper, we will assume that 
the functions m and g are given as in (ig . and use an underlying Clenshaw-Curtis sparse grid. 
For simplicity, we will also only consider isotropic collocation methods, i.e. ai = a 2 = ■ ■ ■ = cxn . 
We then lighten the notation by defining 

Construction of the approximation := requires evaluation of v on a set of collocation 

points Hl C T with cardinality M^. In our case, since the one-dimensional point sets are nested, 
i.e., d' C d'+i for / G N+, so that the multi-dimensional point set used by is given by 


'Hl= [j , 

9(l)=i 

and the nested structure is preserved, i.e., T-Ll C T-Ll+i, to satisfy assumption 13.11 Define the 
difference of the sets A'Hi :=and the number of new collocation points AMl = #{A'Ht,). 
With this nestedness condition, the approximation Xl [u] is a Lagrange interpolating polynomial [3l| , 
and thus dm can be rewritten as a linear combination of Lagrange basis functions. 


Ml 

'^L[v\{y) = '^v{yL,3)'i’LAy) 

i=l 

Ml N 

IG iG{0,l}" n=l 

' -V-' 

’i'L.ji.y) 


(4.5) 


where the index set J{L,j) is defined by 




1 G 


N 

g{l) < L and yL,j G 

n—1 


with i G {0,1} 



For a given L and j, this represents the subset of multi-indices corresponding to the tensor-product 
operators ^ ... ^ c^rn(iN-iN) jjj (|4.2D with the supporting point yLj- Then for each 

1 S J{L,j) and i G {0,1}^, the function Il^Li^ {!,...,m(l„ -*„)},n = 
1,. .., iV, represents the unique Lagrange basis function for the operator 0 . . .^c^rn{iN-iN) 

corresponding to yLj- Therefore, the functions are given by a linear combination of 

tensorized Lagrange polynomials satisfying the “delta property”, i.e., 'itLj'{yL,j) = Sjj' for j,j' = 
1,..., Mi, and is in the required form of (13.31) . 

Finally, to construct the fully-discrete approximation in the space Wh{D) ®Vtp^<s{T) we apply 
the interpolation operator given by (14.5|) . to the semi-discrete solution Uh{x,y) in (13.11) to 

obtain:. 


Ml /Ml \ 

Uh,L{x,y) = XL[uh]{x,y) l^CLJ,iip^ix)\ TLj(y). 

j=l \i=l ) 


(4.6) 
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Due to the delta property of the basis function 'l'ij(y), the interpolation matrix for IlIuh] is a 
diagonal matrix, and thus the coefficient vectors = (cl jp,..., CL,j,Mh) for j = 1 ,..., Ml can 
be computed by independently solving Ml systems of type (EH). 

4.1. Error estimates for fixed L. In what follows, we focus on the linear elliptic problem 
(12.21) described in Examples 12.II and 13.21 and present a detailed convergence and complexity analysis 
of a fully discrete SC approximation, denoted Uh,L, for any fixed level, 1 < E < Emax- As specified 
in Remark 14.11 in this section we consider only the isotropic Smolyak version of SC interpolant given 
by (14.21) . defined on Clenshaw-Curtis abscissas, for solving the parameterized linear elliptic PDE. 
However, our analysis can be extended without any essential difficulty to anisotropic SC methods 
and more complicated underlying PDEs. 

The parameterized elliptic PDE (12.21) admits a weak form that is a symmetric, uniformly coercive 
and continuous bilinear operator on Hq{D); i.e., there exist a,fd>0, depending on a„iin and Omax 
but independent of y, such that for every v,w € Hq{D), 


a{y)\/v Vw dx 


>D 


< «ll«lli/i(D) Ikllifi(D) and / a{y)\Vv\'^ dx. 

D 


In this case, the bilinear form induces a norm ||u||^ = a{y)\Wv\'^ dx, which for functions v{x) = 

c.i4>i{x) G Wh{D), with c = (ci,..., CM;,), coincides with the discrete norm ||c||^(-y), where the 
matrix A{y) is defined in (13.61) . Thus we have 


Continuity: \\c\\j^^y^ = ||u|| < Va , and, (4.7a) 

Ellipticity: y/p ||u||^i(£,) < ||u|| = ||c||^(^) . (4.7b) 

We next state some regularity conditions on the parameterized solution u : T —>■ Hq{D) to the 
parameterized elliptic PDE in Examples 12.II and 13.21 

Assumption 4.1. (Polyellipse analyticity). Let •y = (71,...,y^;) G (l,oo)^, and assume that 
M : P —>■ Hq{D) admits a complex extension u* : —>■ Hq{D), which is analytic on the poly ellipse 


E(7)= n S(n;7n)cC'^, 

l<n<Ar 


where E(n; 7 „) denotes the region bounded by the Bernstein ellipse, 


E (n, 7n) — '^2 A ^ C, I Zji I ^ 7n ^ 


The set £( 7 ) C is the product of ellipses in the complex plane, with foci = ±1, which 
are the endpoints of the domain Pn, n = 1,..., iV. Such ellipses are common in proving convergence 
results for global interpolation schemes. Conditions under which u satisfies Assumption 14.11 can be 
found in [o Theorem 1.2] and [3 Theorem 2.5]. 

In order to investigate the complexity of the fully discrete approximation Uh.L, L G N+, we first 
need to derive sufficient conditions for the error ||m — uu^lWlI to achieve a tolerance of £ > 0, where 
Lg := Lg(r; i7g (D)). Using the triangle inequality, the total error can be split into three parts, i.e., 

||m — Uh.lWl^ ^ Ik ~ UIiWl'^ + \Wh — Uh^bWLl + \Wh,L — uh.lWli- (4-8) 


ei 


es 


14 


62 






The contributions of ei and 62 correspond to the FEM and SC errors, respectively, and have been 
previously examined 31|. The error 63 contributed by the linear solver is often omitted from the 
analysis in the literature, and in practice can be controlled by setting a tight tolerance on the 
iterative solver. However, the analysis presented here is focused on providing cost estimates for the 
iterative solver and reqmres careful consideration of this term. First, we recall error estimates for 
ei and 62 , given from 311. 

Lemma 4.1. Let 7h be a uniform finite element mesh over D C = 1,2,3, with = 

0{\/h'^) grid points. For the random elliptic PDF in Example when u{x, y) € Tg(r; Hq{D) fl 

i4®+^(Zl)),s G N+, the error of the finite element approximation Uh is bounded by 


IIm - Uh\\L 2 < Cfem h®, (4.9) 

e 

where the constant Cfem is independent of h and y. 

Lemma 4.2. Let u satisfy Assumption \4-t\ For L G N“'", the interpolation error u — Xl[u] of 
the sparse grid SC method using Clenshaw-Curtis abscissas can be bounded as 

11^ ~ , (4-10) 

where, for a constant 0 < i5 < 1 , the rate r = (1 — 5) mini<„<Ar log 7 „, and the constant Csc > 0 
depends on N, u, and 6. We remark that the projection of u into the finite element subspace, 
denoted Uh, also satisfies Assumption 14.11 with the same region of analyticity, and therefore the 
application of the interpolant, II, to the semidiscete solution Uh will converge as in ()4.10|1 . 

We now turn our attention to the global solver error 63 in (|4.8I) , which is the error incurred from 
approximating the solution to (1X71) at each sample point. The difference uu^l — Uh,L can be written 
as an interpolant of the solver error, i.e.. 


Uh,L — Uh,L =XL[uh — Uh], 

which represents the solver error amplified by the interpolation operator. For the operator Il[’] in 
(B31) . we have 


< Cl^^ maj^ \\uh{yL,j) '^h{yL,j)\\i^i(^jjy 
Thus, from the ellipticity condition in (I4.7bl) . 


63 < Cl max \\uhiyL,j) - Uh{yL,j, 






< Cl —l= niax | 


CL,j - Clj\ 


Mvl, 






Cl, 


where r is defined to be the tolerance of the linear solver. Note that the expression Uh — Uh is 
only defined at collocation points. The solver error for each fixed yL,j G "Hl is controlled by 
the CG convergence estimate (13.81) . The Lebesgue constant of the operator Il[-] is defined by 
Cl = inaxygrX)^^! L.j{y)\ where 'I'l,^ is given in (14.5|) . We now provide an upper bound of Cl 
in the following lemma. 

Lemma 4.3. The Lebesgue constant for the isotropic sparse-grid interpolation operator Il[-] in 
(lT5l) using the Clenshaw-Curtis rule on T = 0^=1 T" = bounded by 


Cl < [{L + 1){L + 2)f , 


(4.11) 


where L and N are the level of the interpolation operator and dimension of the parameter space, 
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respectively. 

Proof. For each n = 1,.. -jJV, recall that the Lebesgue constants Xi^ of the one-dimensional 
operators are given by 


m{ln) 


i=i 


For Lagrange interpolants based on Clenshaw-Curtis abscissas dUl), we have [l^ 

2 

< - log (m {In) - 1) -I- 1 for In > 2. 


Combining this with the growth rate m{ln) = 2*" ^ + 1 for In > ‘2 given by (14.31) . it is easy to obtain 
that 


Xi„ < 2ln — 1 for In > 2. 

For V S the difference operator for /„ = 1 satisfies 



_ 






< Ai max \v{yn)\. 


For In ^ 2, the triangle inequality yields 




i“(r„) 






< max |z;( 2 /„)|. 


Finally, for v G C'°(r), we bound the interpolant Il[v\ by 


\\Mv]\\ 


L~(r) 


9(l)<i 


L~(r) 


N 


s 12 " E n '• 

g(l)<Ln=l 


max \v 


(L+1 \ ^ 

{yL,j)\ < 2^ ( ^ l\ ^.^max^^ HyL,j)\ 


= [(L + l)(L + 2)]^ max \v{yL,j)l 


which gives the desired estimate. □ 

4.2. Complexity analysis. Now we analyze the cost of constructing Uh.Lma.^, Lmax G N+, with 
the prescribed accuracy e. Here we assume £ > 0 is sufficiently small, and study the asymptotic 
growth of the total costs (13.101) for the accelerated construction of described in ^ For 

comparison, we will also analyze the cost dUD associated with the standard SC method, where 
iterative solvers for the sequence of solutions to the linear systems (1^:71) are seeded with the zero 
vector as an initial guess. According to the error estimates discussed in M.ll a sufficient condition 
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to ensure ||u - < e is that 


||ei|U2 < 

£ 

“ 3’ 

(4.12a) 

\\e 2 \\Ll ^ 

^ *-^SC 6 _ ^ ; 

(4.12b) 

lleallis < 

<(L_+2r^<£ 

(4.12c) 


In section ^we defined K^em and it'acc as the total number of solver iterations used by the standard 
and accelerated SC methods, respectively, to solve (13.711 at each sample point. Now let it'zero(e) 
and itracc(e) represent the minimum values of K^ero and itlacc, respectively, needed to satisfy the 
inequalities (14.121) . Here we aim to estimate upper bounds of KzeTo{£) and Kacds)- Note that, for 
fixed dimension N, level Lmax, and mesh size h, the total number of iterations is determined by the 
inequality ()4.12c|) . Larger values of Lmax and 1/h, lead to higher costs. Thus, the estimation of 
Kzerois) and Kg^ccid) has two steps: (i) Given N and e, estimate the maximum possible h to satisfy 
(I4.12al) and the minimum Lmax that achieves (I4.12bl) ; (ii) Substitute the obtained values into (I4.12cl) 
to estimate upper bounds on iLzero(e) and K^ccd) according to the CG error estimate (13.81) . For 
(i), we have the following lemma, that follows immediately from Lemmas 14.11 and 14.21 

Lemma 4.4. Given the assumptions of Lemmas d.l\ and 14 the error bounds (I4.12ap and 
(l4T^ can be achieved by choosing finite element mesh size h and the sparse-grid level Lmax ac¬ 
cording to 


h{e) 



and Ljjiax(^) 


■ N 

log 2 


log 




(4.13) 


For convenience, we treat the integer quantities K^eroid^ ^"acc(£^), and Lmax(e) as positive real 
numbers in the rest of this section. Now, based on the estimate in Lemma 14.31 for the Lebesgue 
constant CL^axi state the following lemma related to the choice of an appropriate tolerance r(e) 
to satisfy the error bounds (I4.12cl) . 

Lemma 4.5. Let e > 0. Given the assumptions of Lemmas ED and \4.2\ a sufficient condition 
to ensure 63 < e/3 is that 


VPs 

3iL^.ds) + 2 ) 2 ^ ■ 


(4.14) 


Moreover, it holds 


4=(L + 2d^T{e) < Csc for L = 0,.. . ,L„,ax(e) - 1, 

vP 

where L^g^is) is the minimum level given in (14.131) . 

Proof. It is easy to see that (14.141) is an immediate result of (I4.12cl) . For L = 0,..., L^^ds) — 1; 
we have 

4=(L + 2)“r(£) < 4=(L.,„(£) + 2f«r(C < | < < C.. 

which completes the proof. □ 

Using the selected h := h{e), Umax := L^^ds), and r := r(e), we now estimate the upper 
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bounds on the number of CG iterations needed to solve a linear system at a point S ^L^ax- 

To proceed, define 

fczero := max and max kLj for L = 1,..., Tmax, 


where k^j is the number of CG iterations required to achieve \\clj — which, in 

general, depends on the choice of initial vector. Note that, in the case = (0,... ,0)^, there is 
no improvement in the iteration count as the level L increases, so fczero does not depend on L. Now 
we give the following estimates on fczero and {fcacdi”!*- 

Lemma 4.6. Under the conditions of Lemmas and^^ for any G if the CG 

method with zero initial vector is used to solve (l3Jt to tolerance r > 0; then fczero can be bounded by 


fczero 





(4.15) 


Here k = supj^gpK(y), with K,{y) the eondition number of the matrix A{y) corresponding to (13.21) . 
Alternatively, if the initial vector is given by the acceleration method as in (13.51) . then fc^^^, can he 
bounded by 


< log 





(4.16) 


for L — 1, . . . , Lniax- 

Proof. Let yxj be an arbitrary point in Hl, 1 < L < Lmax- Given an initial guess the 
minimum number of GG iterations needed to achieve tolerance r > 0 can be obtained immediately 
from (13.8L that is. 


kL,o = 


log 




•\/^L,j ~ 1 / 


where = A^y^j) is the FE system matrix corresponding to parameter yxj, and klj = i^(yL,j) 
is the condition number of A^j (See Example 13.21) . In the case that = (0,..., 0)^, the estimate 
in (14.151) can be obtained from (I4.7al) . i.e., 


clj - c 


(0) 

'L,j 


— ||cL,il|Ai,3 < V^\\uh\\L°°(r;H^{D)) 


Alternatively, when using Uh,L-i for L = l,...Li„ax to provide initial vectors for the GG solver 
(based on (13.51) 1. for y^j G AHl we use Lemma [43] and (I4.7al) to get the following estimate: 


(0) 

CL.t - (^L,j 


^ (id/i ~ ^^A-l|lL~(r;ffi(D)) fo \Wh,L-l — ^?iT-l|lLoo(r;ffi(D))) 


< -dal Csc e 




n2A, 






< 2^/aCsc 
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This leads directly to the estimate in (14.161) . 0 

In the accelerated case, the sparse-grid interpolant [m/i] must be constructed in the following 
fashion: before solving the system (1X71) corresponding to a sample point yL,j G AT-Ll, we must first 
solve the systems for all sample points in TLl-i- With a total number AMl = #(A'Hl) of new 
linear systems at level L, the total number of CG iterations for the newly added points at level L can 
be bounded by AM^fezero and AM^k^^c^ the standard and the accelerated cases, respectively. 
Then since we find that the total number of iterations for the standard and 

accelerated schemes can be bounded as 

^max 

K^eroie) < hero, and K^cc{e) < ^ AMl k^^c- 

L=1 


This leads to the following estimates. 

Theorem 4.7. Given Assumption E3 and the conditions of Lemmas \f.l\ and \4.iA for e > 0, 
the minimum total number of CG iterations to achieve ||it — Uh,Ln,axl|i| < using zero 

initial vectors is bounded by 


K^eroie) < Cl 


log 


3a, 


N 


C 2 


log 2 


log log 


3a 


iV-l 




log 




a + 2N log log 


1 , /SC,, 


(4.17) 


where k is as defined in Lemma \4.6] and the constants Ci, C 2 , C 3 and C 4 are defined by 

N-l / ^ \ N 

/ p \ 

Gi = 


log 2 


2 

Vn 




*^3 “ ~ (log2) ■ 


(4.18) 


Proof. To achieve the prescribed error, we balance the three error sources that contribute to the 
total error (14.81) . To control ei and 62 , set h = h{e) and Lmax = according to Lemma IXTl 

For the solver error 63 , we choose the solver tolerance r = T(e) according to Lemma [4.51 Then, the 
total number of iterations K^erois) can be bounded by 


Ml„ 


Arzero(e) = X] - 
i=i 


(4.19) 


From Lemma 14.51 and 14.61 we have 




< log 


61/a (Lmax + 2) 

v^e 


2 Af' 


log 


•\/k + 1 
— 1 


(4.20) 
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< 


log 






+ 2N log (Linax + 2) 




^ ^ log ( “ ) + C '4 + 2 A^loglog ( — log 


1 


sa, 


log 


y/n + 1 
^/K — 1 


In addition, following 3l|, Lemma 3.9], we bound the number of interpolation points: 


Ml < V 2^ 

^max — / j 


L=1 


N-l + L 


N - 




< e 


A^—1 nZ/max + 1 


< 2 e 


N-l 


log 


1 + 

sa 




L=1 
N-l 


N-l 

N 


C 2 


log 2 


iV- 1 


log log 


N-l 


^N-1 


(4.21) 


3a, 


N-l 


where in the last line we have used (14.131) to replace Lmax- Substituting (I4.20p and (14.211) into (14.191) 
concludes the proof. □ 

Theorem 4.8. Given Assumption E3 and the conditions of Lemmas ED and \ 4 .‘^ for e > 0, 
the minimum total number of CG iterations it"acc(£); to achieve ||m — u/t,Ln,axll ra < tn Algorithm 
1 , is bounded by 


i^acc(e) < a 


log 


3a, 


1 N 


a + ^^ log log 
log 2 


3a, 


n N-l 


jo + 2 (2 i - l) log (2|i]) + 2JV log log 




3a, 


(4.22) 


where k = supj^gp(K(y)) as in Lemma \4G^ Ci and a are defined as in (|4.18l) . and a is defined by 


a = 2iviog 


2N 
log 2 


a 


+ log(4A/^ ) ■ 


Proof. To achieve the prescribed error, we again choose h = /i(e), Tmax = l^max(£) and t = T(e) 
as in Lemmas 14.41 and 14.51 Then, the total number of iterations i^acc(e) can be bounded by 

Lmax Lmax 

K^ccie)= 

L—1 yL.j^^'Hh L—1 

From Lemma [4.51 and HlBl for L = 1,..., Lujax, we have 


tL < log 


< 


1 
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1 


Hence, 


< 


log 1 

t Vk- 1 / 


1 

log 1 

y/R—l) 


1 

log 1 

fv^+iA 
\/s —1 / 


log 


-in 


L/N ' 


4^ / ^L-Lmaxe 


log (“l/fa...) + rN (2O/O' - 2('-«/») 


log {^\/oi/liCL^A JY / T IN n ^MN\ 

ifacc(e) < ^ ^ I] (2^—- 2(^-1)/'^), 


‘“S (^) 


>‘>8(^) ' 


L = 1 


where S can be bounded using results from geometric sums, i.e., 

^ 2 ^ (2^n,ax/iV _ 2(i-l)/A^) 


L = 1 


< 1 


= e^-i 1 

(- 

< 1 



AT-l Ln 


giV-l j ^ |'2^max/iV _ 2(i-l)/iV^ 2^ 

2-^niax~t”l 2-^max/-^ | 


L =1 


TV- 1 


N-l 


1 - 


1 


2l + l/Af 


2 i+i/Af _ X 


- 2^+-^ max /N 


/ T \ ^ 

M _l_ j ^ 2 ^/^ - 2 ^”“+^ 2 '^”“/^ 


Combining the last two inequalities, along with (I4.2ip . we get 


TVacc(e) < ^ 



s AT-l 

\ oimax + l 

TV- ly 


X 



+ 2Af log (Lmax + 2) + 2rTV ( 2 ^/^ - l) 2^”’“/^ 


Substituting (I4.13P for Lmax concludes the proof. □ 

In the case of the accelerated SC method, an interpolant lL-i[uh], defined by (I4.5p and (I3.4|) . 
must be evaluated for each of the AMl collocation points in AT-Ll- Each interpolant evaluation 
costs 2 TVfi_i — 1 operations, i.e., additions and multiplications, and must be evaluated for each of 
the Mh components of the FEM coefficient vector. Then the interpolation cost on each level is 
TW/jATHj;,(2THl_i — 1) for L = 1,... ,i 2 max(e)- Now we give an estimate of the total interpolation 
cost Cint(e) for our algorithm to achieve the prescribed accuracy e. 

Theorem 4.9. Given Assumption \4-.l\ and the conditions of Lemma EH for sufficiently small 
e > 0 , the total cost of interpolation when using the sparse grid interpolation method in (ITSI) is 


21 































bounded by 


Cint(e)<M,C8(^l0g(^ 


2N 




2(Ar-l) 


where C 2 are defined as in Theorem and Cs = 

Proof. The total interpolation cost is bounded by 

^int {e)<2Mh Y. 

L =2 

<2M, Y 2" )Y 

L=2 ^ ^ Z=1 

-^max(£^} 

'>^1 

iV- 1 


<2M, i: 2u"-i+^)y:2‘ 


L =2 

Lmax{^^ 


Z=1 


<2Mh Y 2^ 


L =2 


A^- 1 + T 
iV- 1 


N-l + l 
A^- 1 

A^ - 1 +; 
N - 1 


■)L-\-l 


<4.Mh 


N 1 + T„iax(e)\ 1 ^L„ax(e) + l 


A^- 1 


< 


16M;,e2(^-i)4^—(") (^1 + 


2(Af-l) 


(4.23) 


Substituting the definition of Tmax(e) from Lemma [4.41 into (14.231) concludes the proof. □ 

Based on Theorems 14.71 HTSl and we finally discuss the savings of the accelerated SC method 

proposed in By comparing the estimates of K^erois) and A"acc(£), we see that the acceleration 
technique reduces log(C' 3 /e) in (14.171) to 2 (2^/^ — l) log(3C'sc/e) in (j4.22l) . Here both terms are 
of the same asymptotic order with respect to e, but the savings from acceleration increases with 
dimension N since (2^/^ — 1) —>■ 0 as A^ —>■ 00. On the other hand, when taking into account the 
cost of interpolation Cint, we must consider the cost Citer of performing each iteration. In the case 
of using CG solvers, Citer is the cost of one matrix-vector multiplication, and will be determined 
by the size of the unknown vector, M/j, and the sparsity of the mass matrix A{y). Thus Citer is 
proportional to the size of the finite element vector, i.e., Citer = where Cd depends on the 

dimension d of the physical domain and choice of finite element basis. For example, without the use 
of a preconditioner, we can assume that the condition numbers of the matrices A{y), for y € T, 
satisfy 


K := sup K{y) < 
yer 

where the constant (7^ > 0 is independent of y £ T Q. Then we can examine the contribution of 
the condition number in Theorems 14.71 and l4. 81 using the inequality log(a;) > {x — l)/x and Lemmas 
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o and 14.41 we bound the terms involving the condition number as 

- 1-^ ^ ^ / 3Cta 

Now as £ —>■ 0, the asymptotic iterative solver costs, Czero = are of the order Mh (y) {log (^) jlog log ( 

while in the accelerated case, the estimate for CdM/j i^acc, is of the same order with respect to £, but 
with an improvement to the constant of — l). For the accelerated method, the additional in¬ 
terpolation costs Cint are of order Mh (log ({■) }^^ (loglog ({■) which is negligible compared 

to the iterative solver complexity. It is clear that, asymptotically, the accelerated method leads to a 
net reduction in computational cost. We remark that for many adaptive interpolation methods, the 
addition of new points already involves evaluation of the current (coarse) interpolant. In this case, 
the cost of interpolation can be ignored, and the accelerated method should be used. 

5. Numerical examples. The goal of this section is to demonstrate the reduction in compu¬ 
tational cost of SC methods using the proposed acceleration technique. In Example 5.1, we first 
use the accelerated SC method to solve an stochastic elliptic PDE with one spatial dimension, and 
compute the overall cost and iteration savings gained by acceleration. Example 5.2 considers a sim¬ 
ilar problem and looks at the number of CG iterations versus the collocation error, comparing the 
implementation of the method using isotropic and anisotropic sparse grids, and demonstrating the 
effect of varying stochastic dimension N on the convergence of the individual systems. In addition, 
as described in Remark 13.31 we extend our acceleration technique to interpolated preconditioners, 
which also exhibit the convergence improvements of the method. Finally, Example 5.3 applies the 
accelerated method to iterative solvers for nonlinear parametrized PDEs. 

The analysis in section 14.11 consisted of two components: (i) estimates for the reduction in 
solver iterations from using acceleration, and (ii) interpolation costs. The interpolation costs can 
be computed exactly for non-adaptive methods, and for adaptive implementations of sparse grid SC 
the interpolation costs can be ignored. In Example 5.1, all error contributions are balanced, and the 
total cost is examined, including both solver iterations and interpolation construction. In Examples 
5.2 and 5.3 we focus only on the number of iterations of the CG solver. 

Example 5.1. We consider the following elliptic stochastic PDE 

-V ■ {a{x,y)Vu{x,y)) = 10 in D x F, , . 

u{x,y) = 0 on dDxT, ' 

where D = [0,1], y = (yi, 2 / 2 , 2 / 3 , r„ = [—1,1], n = 1,..., 4, and the coefficient a is given by: 

log {a{x,y) — 1) = e“^/® (2/1 cos ttx -|- 2/2 sin irx + 2/3 cos 27ra; -|- 2/4 sin 27ra;). (5.2) 

The random variables are independent and identically distributed uniform random variables 

in [—1,1]. In the one-dimensional physical domain, a finite element discretization using linear ele¬ 
ments yields tridiagonal, symmetric positive-definite systems. While this type of system could be 
solved efficiently by direct methods, nevertheless we use CG solvers to demonstrate the convergence 
properties of the acceleration method. 

Table 15.11 compares the standard and the accelerated SG methods, where the error for each 
approximate solution, is computed against a highly refined approximate reference solution 

Uh*,L* with h* = = 10. In Figure [5T] we plot the savings of the accelerated SC method, 

computed according to the cost metrics (13.911 and (13.101) . Since the constants Cfem and Cgc in Lemma 
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HU are not known a priori, to balance the error contributions in (14.121) we use trial and error to 
determine sufficient values h, Lmax, and t to achieve the desired overall error e in the norm. 
Especially for the larger systems, i.e., those with a large number of spatial degrees of freedom, 
significant savings are achieved. The percent savings in the number of iterations versus the cost of 
interpolation are calculated according to 


r — r 


(.^zero ^a.cc') f'int 

.Tf/j Co I^zero 


where Co = 5, since the matrices are tridiagonal. 


Tot. Err 

FE DoFs 

SC Pts 

CG tol 

-^zero 

K 

-*»-acc 

Savings 

1 X 10-2 

255 

137 

1 X 10-3 

28,259 

21,123 

19.4 % 

5 X 10-3 

511 

401 

5 X 10-3 

173,671 

83,884 

42.4% 

1 X 10-3 

2,047 

1,105 

1 X 10-4 

2,001,905 

626,215 

62.3% 

5 X 10-4 

4,095 

2,929 

5 X 10-3 

10,878,352 

1,842,703 

74.5% 

1 X 10-4 

16,383 

7,537 

1 X 10-3 

114,570,175 

12,345,968 

75.1% 


Table 5.1: Comparison in computational cost between the standard and the accelerated SC methods 
for solving (I5.1I) - (I5.2I) . 



Fig. 5.1: Cost (left axis) and percent savings (right axis) of the accelerated SC method versus the 
standard SC method for solving (I5.1I) - (I5.2|) . Costs are computed according to (13.91) and (13.101) . 


Example 5.2. We consider the following stochastic linear elliptic problem 


-V- (a {x,y)Vu {x,y)) 


u{x,y) = 0 
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in D X r, 
on dD X r. 


(5.3) 


















where D = [0,1] x [0,1], r„ = [—\/3, \/3] , n = 1,..., and x = {xi,X 2 ) is the spatial variable. The 
random diffusion term has one-dimensional spatial dependence given by 


log(a(x, y) - 0.5) = 1 + yi 4- ^ CnV^n(^)yn, 


N 


n—2 


(5.4a) 


where 


c„ := (^i?)i/ 2 exp (Ln/2j7ri?)^ ^ ^ 


n > 1 


(5.4b) 


and 


‘Pnix) := 


, f [n/2j TTXi \ 
Rp )' 
I n/2 1 TTXi \ 

' IL ) 


n even 


n odd. 


(5.4c) 


The random variables {yn}n=i i.i.d. and are each uniformly distributed in [—-\/3, v^], with zero 
mean and unit variance, i.e., E[y„] = 0, and £[?/„?/„] = Snm, for n,m € N+. The finite dimensional 
stochastic diffusion a represents the A^-term truncation of an expansion of a random field with 
stationary covariance function, given by 


Cov [log (a — 0.5)] {xi,x'i) = exp 


{xi - x\f 

Ri 


(5.5) 


where xi^x'i G [0,1], and Rc is the physical correlation length for the random field a. The parameter 
Rp in (I5.4cl) is given by Rp = max{l, 2Rc} and R is given by i? = RdRp. Then and (pn(x) are the 
eigenvalues and eigenfunctions associated with (15.5|) . Here we will consider two correlation lengths, 
namely Rc = 1/2, and Rc = 1/64, where Figure [Ol shows the corresponding decay of eigenvalues. 
For the spatial discretization, we use a finite element approximation on a regular triangular mesh 
with linear finite elements and 4225 degrees of freedom. The CG method is used for the linear solver 
with diagonal preconditioners and a tolerance of 10“^^. 

First, for Rc = 1/64, the error and total iteration count of both the standard case, using 
zero initial vectors, and accelerated SC construction, computed using several dimensions N, are 
summarized in Table [521 The error is measured using the expectation of the approximate solutions, 
- E[m/i,l*]|1l2 (£,), for Lmax = 1, • ■ •, 7, where the “exact” solution ¥\uh,LA is computed 
using L* = 8. We compare these errors against the cumulative total number of iterations, ATzero and 
ATacc, needed to construct 

An alternative approach to accelerating SC methods is found in 2l|. For a particular SC 
level Lmax, this method orders the collocation points lexicographically, with each dimension ordered 
according to the decay of the eigenvalues in (I5.4al) . We also implemented a similar method without 
the sequential ordering; for a given level L, at each new collocation point in ATL^ the solution at the 
nearest collocation point from lower levels is given as an initial guess to accelerate the CG solver. 
We refer to this method as the “nearest neighbor” approach. Figure [531 shows the average number 
of iterations needed to solve the linear system (EH), where the average is taken over the new points 
at level L, i.e., AR^, for L = 1,... ,7. We compare our interpolated acceleration algorithm, the 
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Fig. 5.2: First 19 eigenvalues for (15.5|) for correlation length Rc = 1/64,1/2. 



Error 

SC Pts 

Kzero 

^acc 

Savings in K 

N=3 

3.83e-8 

25 

6,780 

5,991 

11.6% 

9.57e-10 

69 

18,893 

14,628 

22.6% 

9.86e-12 

177 

48,691 

27,765 

43.0% 

N=5 

5.28e-07 

61 

17058 

15095 

11.6% 

1.03e-08 

241 

67,955 

53,992 

20.6% 

1.44e-10 

801 

226,597 

150,241 

33.7% 

N=7 

2.43e-08 

589 

168,237 

136,072 

19.1% 

6.63e-10 

2,465 

706,049 

500,718 

29.1% 

1.94e-ll 

9,017 

2,585,970 

1,496,391 

42.1% 

N=9 

1.68e-07 

1,177 

338,428 

277,583 

18.0% 

7.83e-09 

6,001 

1,729,337 

1,273,895 

26.3% 

8.86e-ll 

26,017 

7,505,343 

4,719,820 

37.1% 

N=ll 

2.59e-07 

2,069 

596,368 

495,705 

16.9% 

2.43e-08 

12,497 

3,608,185 

2,736,615 

24.2% 

1.95e-09 

63,097 

18,231,420 

12,139,658 

33.4% 


Table 5.2: Iteration counts and savings of the accelerated SC method for solving (j5.3ll - ()5.4ll with 
correlation length Rc = 1/64, and stochastic dimensions iV = 5, 7,9, and 11. 


nearest neighbor approach, and standard SC method without acceleration, for N = 3 and = 11, 
using Rc = 1/64. The interpolated initial vector provided by the acceleration algorithm yields a 
reduction in the average number of iterations at each level, which increases with L. Figure [531 also 
shows the effect of using the nearest neighbor solution as the initial vector, which provides some 
improvement over the standard case using zero initial vectors, but the savings do not match those of 
our acceleration scheme. Note that since the number of new collocation points grows exponentially 
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Fig. 5.3: Comparison of the average CG iterations per level for solving problem (I5.3D - (I5.4D with 
dimensions TV = 3 (left) and = 11 (right), and correlatioir length Rc = 1/64. 


with each level (cf (14.31) 1. there is an increase in total iteration savings over successive levels in both 
the nearest neighbor and accelerated case. 

The left plot of Figure 15.41 shows the total iteration savings achieved by the acceleration algo¬ 
rithm with different maximum collocation levels Lmax = 1,..., 6. The savings are measured as the 
percentage reduction in the cumulative iteration count up to level Lmax, relative to standard case 
using zero initial vectors, i.e., (iVzero — Kzero)/Kzero- Here we also see the effect of stochastic dimen¬ 
sion on the convergence of SC methods: as N increases, our algorithm provides less accurate initial 
guesses for a given maximum SC level Lmax- This can also be seen by comparing the left and right 
plots of Figure 15751 which show how the average number of iterations at a given SC level L changes 
from TV = 3 to TV = 11. On the other hand, the right plot of Figure [574] shows the same total iteration 
savings now plotted versus error. As above, the error is measured as — ^[uh,L*]\\L^(D)j 

with L* = 7. These results are in agreement with the theoretical asymptotic estimates from Theorem 
14.81 which predict an increased savings vs error for larger dimensions. 

Next we examine the effect of the correlation length, i?c, on our acceleration algorithm. Larger 
correlation lengths result in faster decay of eigenvalues of the covariance function (15.51) (see Figure 
15.2p . and implies that u{y) depends on certain components of the vector y more than others, which 
reduces the effectiveness of isotropic methods. Figure 15.51 plots the convergence of the error in 
versus the total number of CG iterations for TV = 3 and TV = 11, and for both Rc = 1/2 
and Rc = 1/64. The larger correlation length, Rc = 1/2, results in slower convergence of the SC 
interpolant than for Rc = 1/64, but note that the accelerated method reduces the total iteration 
count in both cases. 

On the other hand, we can employ anisotropic methods to increase the efficiency of SC in the 
case of larger correlation lengths [30|. Anisotropic SC methods will place more points in directions 
corresponding to large eigenvalues of (15.51) . and the importance of each dimension is encoded in a 
weight vector (see (I4.3|) b Figure 15751 plots the average number of iterations for problem (I5.3I) - (I5.4I) 
with a relatively large correlation length R^ = 1/2, and TV = 11. Here we employ the weights given 
by an a posteriori selection described in [30|, i.e., the weight vector a G K.'^, with ai = 0.85, a 2 = 
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Fig. 5.4: Percentage cumulative reduction in CG iterations vs level (left) and error (right) for 
solving (I5.3p - (l5.4p using our accelerated approach, with N = 5,7,9, and 11 and for correlation 
length Rc = 1/64. 




10 


10 10 
Total CG Iterations 


10 


Fig. 5.5: The convergence of the SC approximation for solving (I01)-(I01). using CG, with and 
without acceleration, for correlation lengths Rc = 1/64,1/2, and dimensions N = 3 (left), and 
iV = 11 (right). 


as = 0 . 8,04 = 05 = 1.0,06 = 07 = 1 . 6,08 = 0:9 = 2.6, oip = on = 3.7. The acceleration method 
decreases the average number of iterations needed to solve the linear system, but the effect is not 
as pronounced as in the case of an isotropic SC method. This occurs because the isotropic method 
places far too many points in relatively unimportant directions, thus the dependence of u{y) on 
a certain component t/„ of y may be well approximated at very low levels. Anisotropic methods 
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Fig. 5.6: Average CG iterations per level for solving problem (I5.3|) - (I5.4|) for = 11 and with 
correlation length Rc = 1/2, using an isotropic SC (left) and anisotropic SC (right). The inefhciencies 
from using an isotropic grid are partially offset by increased gains from acceleration. 


exhibit better convergence with respect to (and lower interpolation costs) versus isotropic 

methods, yet we see here that the acceleration algorithm helps to somewhat offset the inefhciency 
of isotropic methods for anisotropic problems. 

In the preceding results we have used a simple diagonal preconditioner strategy. As described in 
Remark 13.31 we can also construct efficient preconditioners with our acceleration scheme. Table [5731 
shows the effectiveness of the preconditioning strategy for solving equations (I5.3|1 - (I5.4|1 . with N = 7 
and Rc = 1/64, where we compare the average number of iterations needed to solve (j3.7l) at each new 
point hlj G A'Hi at a given level L. Here we compute an incomplete Cholesky preconditioner for 
each linear system on the levels L = 1,..., Tpc, for Lpc = 1, 2, and 3, and use these to provide an 
“accelerated” preconditioner (13.111) for the systems on the remaining levels Lpc + 1, • ■ ■, Tmax- We 
compare this against the cases where a simple diagonal preconditioner and an incomplete Cholesky 
preconditioner are used for each system. The three-level accelerated preconditioner reduces the 
average number of iterations to within a decimal point of the incomplete Cholesky preconditioner, 
and the cost of computing the low-level preconditioners and interpolating is relatively cheap in 
comparison. 

Example 5.3. The preceding experiments demonstrate the benefits of using acceleration to 
improve the convergence of individual iterative linear solvers. In the case of a nonlinear PDE, the 
possibilities for savings can be even greater than the linear cases above, since convergence of a 
nonlinear solver may be slow or even unattainable from a poor initial vector. In this example, we 
consider the problem 


{ -V • (a {x, y) Vu {x, y)) + F[u]{x, y) = x in H x F, 

u(0,y) = 0 in r, 

u'{l,y) = l in r. 


where a is given by (Ol) . D = [0,1], r„ = [—1,1], n = 1,..., 4, and F'[it] is some nonlinear function 
of u. In what follows, we consider the nonlinear functions E[m] = rt®, and F'[u] = uu'. 
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GG iterations for standard SG 

Level 

No PC 

Diag PC 

Inc. Choi. 

Lpc = 1 

o 

II 

to 

Lpc = 3 

1 

243 

243 

55 

55 

- 

- 

2 

311.8 

278.4 

54.7 

60.7 

54.7 

- 

3 

332.3 

284.9 

54.6 

63.5 

54.9 

54.6 

4 

341.0 

286.1 

54.6 

65.2 

55.3 

54.6 

5 

345.8 

286.7 

54.6 

66.2 

55.5 

54.6 

6 

348.4 

286.9 

54.6 

66.7 

55.6 

54.6 

GG iterations for accelerated SC 

Level 

No PC 

Diag PC 

Inc. Choi. 

Lpc = 1 

Lpc = 2 

Lpc = 3 

1 

243 

243 

55 

55 

- 

- 

2 

299.3 

264.6 

52.9 

58.4 

52.9 

- 

3 

295.8 

251.3 

49.1 

57.1 

49.4 

49.1 

4 

270.8 

225.8 

43.7 

52.3 

44.2 

43.7 

5 

237.0 

194.3 

37.3 

45.8 

38.0 

37.3 

6 

186.1 

151.9 

28.9 

36.0 

29.5 

28.9 


Table 5.3: Average iteration counts for the standard SC method (top), and the accelerated SC 
method (bottom) using six preconditioner schemes to solve (I5.3I) - (I5.4I) with N = 7, and Rc = 1/64. 
From left to right: no preconditioner, diagonal preconditioners, incomplete Cholesky preconditioners, 
and accelerated preconditioners (13.111) built using incomplete Cholesky preconditioners with Lpc = 
1,2,3. 


Nonlinear problems are typically solved with the use of iterative methods such as Picard itera¬ 
tions or Newton’s method. We implement a combination of these methods that begins with Picard 
iterations, then utilizes Newton’s method once the relative errors are small. For spatial discretiza¬ 
tion, we use piecewise linear finite elements on [0,1] with a mesh size oi h = 1/500, and solved the 
resulting systems at each iteration using exact methods. The stopping criterion for the solver is a 
relative tolerance of 10“® in the P norm. 

Results for these experiments are given in Figure [5?7l For each SC level, L = 1,..., 8, we plot 
the average number of nonlinear iterations, where the average is taken over the set of points which 
are new to level L, namely M-Ll- Finally, we show the total computational time in Table 1?^ 
for different maximum levels of stochastic approximation, measured on a workstation with 1.7GHz 
dual core processors and 8 GB of RAM. We note that in Table 15.41 the size of the finite element 
system is fixed. Thus, as we move to higher levels of collocation, the stochastic approximation 
becomes relatively more expensive to compute compared to the solving the finite element systems. 
This is why the savings begin to decrease after level 5, even though Figure 15.71 shows dramatic 
savings in iterations for higher levels. Furthermore, the reason for the negative savings for a level 
L = 2 stochastic approximation is that the interpolant is not yet accurate enough to overcome the 
additional cost of the acceleration. 

6. Conclusion. In this work, we proposed and analyzed an acceleration method for construc¬ 
tion of sparse interpolation-based approximate solutions to PDFs with random input parameters. 
The acceleration method exploits the sequence of increasingly accurate approximate solutions to 
provide increasingly good initial guesses for the underlying iterative solvers that are used at each 
sample point. We have developed this method using a global Lagrange polynomial basis but the 
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Fig. 5.7: Average number of Newton iterations per level for solving problem (|5.6p with F[u] = uu' 
(left) and F[u\ = (right). 


SC Level 

2 

3 

4 

5 

6 

F[m] = acc 

.03018 

.113832 

.2746 

.7039 

2.33314 

Ejit] = u®, zero 

.025976 

.119256 

.339678 

.949184 

2.61958 

% Savings 

-16.2 

4.5 

19.2 

25.8 

10.9 

F[u] = uu', acc 

.027754 

.089082 

.22706 

.629451 

2.05741 

F[u] = uu', zero 

.026527 

.090435 

.273355 

.895027 

2.4008 

% Savings 

-4.6 

1.5 

16.9 

29.7 

14.3 


Table 5.4: Computational time in seconds for computing solution to problem 15.61 


method can easily be extended to other non-intrusive methods. 

While our method takes advantage of the natural structure provided by hierarchical SC methods, 
we do not take advantage of any hierarchy in the spatial approximation. As explained in Remark l3.ll 
our method may be used in combination with the multilevel method to accelerate the construction of 
stochastic operators, and reuse information from level to level. The combination of the acceleration 
scheme with multilevel methods will be the subject of future work. 

We rigorously studied error estimates in the special the case of linear elliptic PDEs with ran¬ 
dom inputs, providing complexity estimates for the proposed method. Several numerical examples 
confirm the expected performance. While the analysis of >14.11 applies to linear stochastic PDEs, 
the acceleration method may be even more well suited to nonlinear problems, as convergence rates 
may be improved, based on the choice of a good initial guess for nonlinear iterative solvers. A final 
numerical example demonstrates the advantage of our approach to nonlinear problems. A more 
rigorous study of acceleration for nonlinear solvers and extension to time dependent problems may 
provide interesting opportunities in the future. 
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